      SUBROUTINE PLOT_FINAL_PDFS( ICSET, IUNIT, JUNIT, AERO )
      USE AERO_PARAM, ONLY: PI6, DENSP, IXXX, IYYY, ILAY
      USE AERO_CONFIG, ONLY: NMODES, NAEROBOX
      USE AERO_SETUP, ONLY: SIG0, CONV_DPAM_TO_DGN, NUMB_MAP, MODE_NAME
      USE AERO_DIAM
      IMPLICIT NONE

      ! Arguments.
 
      INTEGER, INTENT(IN) :: ICSET         ! identifies test case [1]
      INTEGER, INTENT(IN) :: IUNIT         ! logical unit number for number pdf [1]
      INTEGER, INTENT(IN) :: JUNIT         ! logical unit number for mass   pdf [1]
      REAL(8), INTENT(IN) :: AERO(NAEROBOX)! aerosol conc. [ug/m^3] or [#/m^3]

      ! Local variables. 

      INTEGER :: I, N
      INTEGER, PARAMETER :: NBINS = 30! 200    ! number of bins [1]      
      REAL(8) :: DGRID(NBINS)              ! fixed diameter        grid [um]
      REAL(8) :: MGRID(NBINS)              ! fixed mass/particle   grid [ug/particle]
      REAL(8) :: DLOWER(NBINS)             ! lower boundary fixed diameter grid [um]
      REAL(8) :: DUPPER(NBINS)             ! upper boundary fixed diameter grid [um]
      REAL(8) :: NTOT(NMODES)              ! number concentration for each mode [#/m^3]
      REAL(8) :: PDF(NBINS,2,NMODES)       ! number or mass conc. at each grid point [#/m^3] or [ug/m^3]       
      REAL(8) :: DNDLOGD(NMODES)           ! dN/dlog10(Dp) [ #/m^3]
      REAL(8) :: DMDLOGD(NMODES)           ! dM/dlog10(Dp) [ug/m^3]
      REAL(8) :: RDMIN                     ! reciprocal of DMIN to optimize coagulation [1/um]
      REAL(8) :: RDLOGDSC                  ! reciprocal of log10 of the grid spacing [1]
      REAL(8) :: SCALE, F, SUM1, SUM2      ! scratch variables 
      REAL(8) :: DMINL, DMAXL, DG          ! diameters [um]  
      REAL(8) :: FLN                       ! function for lognormal distribution [1]  
      REAL(8), PARAMETER :: DMIN =  0.001D+00   ! smallest particle diameter of the discrete grid [um]
      REAL(8), PARAMETER :: DMAX = 20.000D+00   ! largest  particle diameter of the discrete grid [um]

      SELECT CASE( ICSET )
      CASE ( 10 )
        DMAXL = 0.600D+00
        DMINL = 0.006D+00
      CASE ( 11 )
        DMAXL = 1.000D+00
        DMINL = 0.001D+00
      CASE ( 12 )
        DMAXL = DMAX
        DMINL = DMIN
      CASE ( 13 )
        DMAXL = DMAX
        DMINL = DMIN
      CASE ( 14 )
        DMAXL = 100.000D+00
        DMINL =   0.001D+00
      CASE ( 15 )
        DMAXL = DMAX
        DMINL = DMIN
      CASE ( 16 )
        DMAXL = DMAX
        DMINL = DMIN
      CASE ( 17 )
        DMAXL = 200.000D+00
        DMINL =   0.001D+00
      CASE ( 18 )
        DMAXL = 8.000D+00
        DMINL = 0.001D+00
      CASE ( 19 )
        DMAXL = DMAX
        DMINL = DMIN
      CASE ( 20 )
        DMAXL = DMAX
        DMINL = DMIN
      CASE DEFAULT
        DMAXL = DMAX
        DMINL = DMIN
        ! WRITE(*,*)'Should not reach CASE DEFAULT in subr. plot_final_pdfs.'
        ! STOP
      END SELECT

      SCALE    = ( DMAXL / DMINL )**(1.0D+00/REAL(NBINS-1))
      RDLOGDSC = 1.0D+00 / LOG10( SCALE )
      RDMIN    = 1.0D+00 / DMINL
      DO I=1, NBINS
        DGRID(I)  = DMINL * SCALE**(I-1)                  ! [um]
        DLOWER(I) = DGRID(I) / SCALE**0.5D+00             ! [um]
        DUPPER(I) = DGRID(I) * SCALE**0.5D+00             ! [um]
        MGRID(I)  = 1.0D-06 * DENSP * PI6 * DGRID(I)**3   ! [ug/particle]
        DO N=1, NMODES
          DG = 1.0D+06 * DIAM(IXXX,IYYY,ILAY,N) * CONV_DPAM_TO_DGN(N)   ! convert [m] to [um] and Dbar to Dg
          NTOT(N) = AERO( NUMB_MAP(N) )
          F = NTOT(N) * FLN( DGRID(I), DG, SIG0(N) )
          PDF(I,1,N) = F * ( DUPPER(I) - DLOWER(I) )
          PDF(I,2,N) = PDF(I,1,N) * MGRID(I)
          DNDLOGD(N) = PDF(I,1,N) * RDLOGDSC * 1.0D-06           ! convert from [#/m^3] to [#/cm^3]
          DNDLOGD(N) = MAX( DNDLOGD(N), 1.0D-30 )
          DMDLOGD(N) = PDF(I,2,N) * RDLOGDSC                     ! [ug/m^3]
          DMDLOGD(N) = MAX( DMDLOGD(N), 1.0D-30 )
          !------------------------------------------------------------------------------------------------------------
          ! WRITE(IUNIT, *) I,N, NUMB_MAP(N), NTOT(N)
          ! WRITE(IUNIT,92) I,N, DGRID(I),DG,SIG0(N),NTOT(N),F,PDF(I,1,N)
          !------------------------------------------------------------------------------------------------------------
        ENDDO
        WRITE(IUNIT,91) I, DGRID(I), DNDLOGD(:)
        WRITE(JUNIT,91) I, DGRID(I), DMDLOGD(:)
      ENDDO

      DO N=1, NMODES
        SUM1 = 0.0D+00
        SUM2 = 0.0D+00
        DO I=1, NBINS
          SUM1 = SUM1 + PDF(I,1,N)
          SUM2 = SUM2 + PDF(I,2,N)
        ENDDO
        ! WRITE(*,'(I4,A6,F12.1,F12.6)') N, MODE_NAME(N), 1.0D-06*SUM1, SUM2
      ENDDO


91    FORMAT(I5,17D13.5)
92    FORMAT(2I5,7D13.5)
      RETURN 
      END SUBROUTINE PLOT_FINAL_PDFS


      SUBROUTINE TEST_GETACTFRAC
      USE AERO_ACTV, ONLY: GETACTFRAC 
      IMPLICIT NONE
      INTEGER, PARAMETER :: NMODES5 = 5  ! number of activating modes [1]      
      REAL(8) :: TKELVIN                 ! absolute temperature [K]
      REAL(8) :: PTOT                    ! ambient pressure [Pa]
      REAL(8) :: WUPDRAFT                ! updraft velocity [m/s]
      REAL(8) :: XNAP5(NMODES5)          ! number concentration of particles for each mode [#/m^3]
      REAL(8) :: XMAP5(NMODES5,5)        ! number concentration of particles for each mode [#/m^3]
      REAL(8) :: RG5(NMODES5)            ! geometric mean radius of the lognormal for each mode [um]
      REAL(8) :: SIGMAG5(NMODES5)        ! geometric standard deviation of the lognormal for each mode [um]
!      REAL(8) :: BIBAR5(NMODES5)         ! hygroscopicity parameter for each mode [1]
      REAL(8) :: FRACACTN5(NMODES5)      ! activating fraction of the aerosol number conc. for each mode [1]
      REAL(8) :: FRACACTM5(NMODES5)      ! activating fraction of the aerosol mass   conc. for each mode [1]
      REAL(8) :: AC5(NMODES5)            ! minimum dry radius for activation for each mode [um]
      REAL(8) :: NACT5(NMODES5)          ! activating number concentration for each mode [#/m^3]
      REAL(8) :: MACT5(NMODES5)          ! activating mass   concentration for each mode [ug/m^3]

      ! Test routine GETACTFRAC. 

      XNAP5(1)   =    200.00D+06  ! [#/m^3] sulfate
      XNAP5(2)   =    100.00D+06  ! [#/m^3] BC
      XNAP5(3)   =    100.00D+06  ! [#/m^3] OC
      XNAP5(4)   =     20.00D+06  ! [#/m^3] dust
      XNAP5(5)   =      1.00D+06  ! [#/m^3] sea salt
      XMAP5(:,:) =      0.00D+00  ! [ug/m^3] 
      XMAP5(1,1) =      1.00D+00  ! [ug/m^3] 
      XMAP5(2,2) =      1.00D+00  ! [ug/m^3] 
      XMAP5(3,3) =      1.00D+00  ! [ug/m^3] 
      XMAP5(4,4) =      1.00D+00  ! [ug/m^3] 
      XMAP5(5,5) =      1.00D+00  ! [ug/m^3] 
      RG5(1)     =      0.08D+00  ! [um] 
      RG5(2)     =      0.08D+00  ! [um] 
      RG5(3)     =      0.08D+00  ! [um] 
      RG5(4)     =      0.30D+00  ! [um] 
      RG5(5)     =      1.00D+00  ! [um] 
      SIGMAG5(:) =      1.80D+00  ! [1] 
      TKELVIN    =    294.00D+00  ! [K]
      PTOT       = 100000.00D+00  ! [Pa]
      WUPDRAFT   =      0.50D+00  ! [m/s]

      CALL GETACTFRAC(NMODES5,XNAP5,XMAP5,RG5,SIGMAG5,TKELVIN,PTOT,WUPDRAFT,AC5,FRACACTN5,FRACACTM5,NACT5,MACT5)

      ! WRITE(*,91) FRACACTN5(:)
      ! WRITE(*,91) FRACACTM5(:)

91    FORMAT(5F12.8)
      RETURN
      END SUBROUTINE TEST_GETACTFRAC


      SUBROUTINE TEST_DEPV
      USE AERO_PARAM, ONLY: IXXX, IYYY, ILAY
      USE AERO_DEPV, ONLY: GETDEP_V, GET_AERO_DEPV, VDDEP_AERO
      IMPLICIT NONE
      INTEGER, PARAMETER :: N = 1
      INTEGER :: I
      REAL(8) :: DBLKTA       ! air temperature [ K ]
      REAL(8) :: DBLKDENS     ! air density  [ kg/m**3 ]
      REAL(8) :: DXLM         ! atmospheric mean free path [ m ]
      REAL(8) :: DAMU         ! atmospheric dynamic viscosity [ kg/(m s) ]
      REAL(8) :: DBLKWSTAR    ! convective velocity scale [ m/s ]
      REAL(8) :: DBLKUSTAR    ! friction velocity [ m/s ]
      REAL(8) :: DBLKRA       ! aerodynamic resistance [ s/m ]
      REAL(8) :: DGN_DDEP(N)  ! 
      REAL(8) :: XLS_DDEP(N)  ! 
      REAL(8) :: DEN_DDEP(N)  ! 

      DBLKTA    = 288.15D+00      ! [ K ]
      DBLKDENS  = 1.225D+00       ! [ kg/m**3 ]
      DXLM      = 6.6332D-08      ! [ m ]
      DAMU      = 1.7894D-05      ! [ kg/(m s) ]
      DBLKWSTAR = 1.0D+00         ! [ m/s ]
      DBLKUSTAR = 0.5D+00         ! [ m/s ]
      DBLKRA    = 5.0D+00         ! [ s/m ]
      DGN_DDEP(:) = 0.001D+00     ! [um]
      XLS_DDEP(:) = LOG(1.6D+00)  ! [1] 
      DEN_DDEP(:) = 1.5D+00       ! [g/cm^3] 
      IXXX = 1
      IYYY = 1

      WRITE(*,'(A/)')' I, DG(um), VDEPn (m/s), VDEPm (m/s) for Sigma = 1.6'
      DO I=1, 26
        CALL GET_AERO_DEPV(N,DBLKTA,DBLKDENS,DXLM,DAMU,DBLKWSTAR,DBLKUSTAR,DBLKRA,DGN_DDEP,XLS_DDEP,DEN_DDEP)
        WRITE(*,91) I, DGN_DDEP(1), VDDEP_AERO(IXXX,IYYY,1,1), VDDEP_AERO(IXXX,IYYY,1,2)
        DGN_DDEP(:) = DGN_DDEP(:) * 10.0**(0.2)
      ENDDO

91    FORMAT(I4,F9.4,2X,2F11.6)
      RETURN
      END SUBROUTINE TEST_DEPV


      REAL(8) FUNCTION FLN(X,XG,SIGMAG)
      REAL(8) :: X      ! particle radius or diameter [any units]
      REAL(8) :: XG     ! geometric mean radius or diameter [any units]
      REAL(8) :: SIGMAG ! geometric standard deviation [monodisperse = 1.0]
      REAL(8), PARAMETER :: SQRTTWOPI = 2.506628275D+00
      FLN = EXP(-0.5D+00*(LOG(X/XG)/LOG(SIGMAG))**2) / (X*LOG(SIGMAG)*SQRTTWOPI)
      RETURN
      END FUNCTION FLN




